Secular Evolution of Hierarchical Planetary Systems 
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ABSTRACT 

We investigate the dynamical evolution of coplanar hierarchical two-planet systems 
where the ratio of the orbital semimajor axes a = ai/a2 is small. Hierarchical two- 
planet systems are likely to be ubiquitous among extrasolar planetary systems. We 
show that the orbital parameters obtained from a multiple Kepler fit to the radial ve- 
locity variations of a host star are best interpreted as Jacobi coordinates and that Jacobi 
coordinates should be used in any analyses of hierarchical planetary systems. An ap- 
^ , proximate theory that can be applied to coplanar hierarchical two-planet systems with a 

I wide range of masses and orbital eccentricities is the octupole-level secular perturbation 

^) ' theory, which is based on an expansion to order and orbit-averaging. It reduces the 

■ coplanar problem to one degree of freedom, with ei (or 62) and zui — W2 as the relevant 

■ phase-space variables (where ei^2 are the orbital eccentricities of the inner and outer 
Q . orbits and wi^2 are the longitudes of periapse). The octupole equations show that if 

r~| . the ratio of the maximum orbital angular momenta, A = L1/L2 ~ (mi/m2)a^''^, for 

given semimajor axes is approximately equal to a critical value Acriti then libration of 
O \ w\ — VJ2 about either 0° or 180° is almost certain, with possibly large amplitude varia- 

' tions of both eccentricities. From a study of the HD 168443 and HD 12661 systems and 

■ their variants using both the octupole theory and direct numerical orbit integrations, 
we establish that the octupole theory is highly accurate for systems with a < 0.1 and 

^ . reasonably accurate even for systems with a as large as 1/3, provided that a is not too 

close to a significant mean-motion commensurability or above the stability boundary. 
The HD 168443 system is not in a secular resonance and its wi — W2 circulates. The 
HD 12661 system is the first extrasolar planetary system found to have wi — W2 librat- 
ing about 180°. The secular resonance means that the lines of apsides of the two orbits 
are on average anti-aligned, although the amplitude of libration of wx — W2 is large. 
The libration of w\ — 'W2 and the large-amplitude variations of both eccentricities in 
the HD 12661 system are consistent with the analytic results on systems with A ~ Acrit- 
The evolution of the HD 12661 system with the best-fit orbital parameters and sini = 1 
(i is the inclination of the orbital plane from the plane of the sky) is affected by the close 
proximity to the 11:2 mean- motion commensurability, but small changes in the orbital 
period of the outer planet within the uncertainty can result in configurations that are 
not affected by mean-motion commensurabilities. The stability of the HD 12661 system 
requires sini > 0.3. 
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1. INTRODUCTION 

Extrasolar planet searches using high-precision radial velocity observations to measure the 
reflex motion of the host stars along the line of sight have now yielded over 100 extrasolar planets.^ 
The discoveries include 10 multiple planet systems with either two or, in the cases of v And and 
55 Cnc, three detected planets. There are also indications that about half of the stars with one 
known planet are likely to have additional detectable distant companions (Fischer et al. 2001). 
Many of the known multiple planet systems exhibit interesting dynamics. The two planets about 
GJ 876 are deep in three orbital resonances at the 2:1 mean- motion commensurability (Laughlin 
k, Chambers 2001; Lee & Peale 2002), and there are possibly two other systems with mean-motion 
resonances: 2:1 for the two planets about HD 82943^ (Gozdziewski & Maciejewski 2001) and 3:1 
for the inner two planets about 55 Cnc (Marcy et al. 2002; Lee &; Peale 2003). The outer two 
planets of the v And system are apparently locked in a secular resonance with the lines of apsides of 
the two orbits being aligned on average (Rivera &: Lissauer 2000; Lissauer & Rivera 2001; Chiang, 
Tabachnik, &; Tremaine 2001), and as we shall demonstrate in this paper, the two planets of the 
HD 12661 system are also locked in a secular resonance, but with the lines of apsides being anti- 
aligned on average. The resonances are often vital in ensuring the long term stability of the multiple 
planet systems. 

Although the origin, evolution, and stability of the orbital configurations of multiple planet 

systems can usually be analyzed with direct numerical integrations of the full equations of motion, 
a theoretical understanding of the dynamics often requires the development and application of 
theories with analytic approximations. The approximate theories also allow one to explore the 
parameter space much more rapidly than direct numerical integrations in the regions where the 
approximations are valid. 

Because the orbital eccentricities and inclinations of the planets in the solar system are gen- 
erally small while the ratios of the orbital semimajor axes of adjacent planets are generally large 
(oj/aj+i > 0.5 except for the Mars-Jupiter pair), the classical perturbation theory developed for 
the solar system is based on an expansion of the disturbing functions that describe the mutual grav- 
itational interactions of the planets in powers of the eccentricities and inclinations. In particular, 
the Laplace-Lagrange secular solution for the evolution of a two-planet system that is not affected 
by mean-motion commensurabilities is based on retaining just the secular terms (i.e., those not 
involving the mean longitudes) in the disturbing functions up to second order in the eccentricities 
and inclinations (see, e.g., Murray & Dermott 1999). Since high orbital eccentricities are common 
among extrasolar planets and the distribution of the orbital eccentricities of the extrasolar planets 
in the known multiple planet systems is similar to that of the single planets (Fischer et al. 2003), 
the classical Laplace-Lagrange secular perturbation theory is not in general adequate for describing 



^See, e.g., http://www.obspm.fr/planets for a continuously updated catalog. 
^See http://obswww.unige.ch/~udry/planet/lid82943syst.litml. 
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the secular evolution of the extrasolar planetary systems that are not affected by mean-motion 
commensurabilities. 

On the other hand, many of the known multiple planet systems are hierarchical in the sense 
that the ratio(s) of the scmimajor axes (01/02 for the two-planet systems and 01/02 and/or 02/03 
for the three-planet systems) are small. There are 12 pairs of adjacent planets in the 10 known 
multiple planet systems (two of which have three planets). Among the 9 pairs that are not known 
or suspected to be in mean-motion resonances, five have Oj/oj+i < 0.1 and all but one (47 UMa) 
have aj/aj+i < 1/3. Thus a secular perturbation theory for a two-planet system that is based on 
an expansion in 01/02 and valid for high eccentricities should provide an accurate description of 
the secular evolution of many extrasolar planetary systems. 

It is useful to consider hierarchical two-planet systems in the context of the general hierarchical 
triple systems in which a third body orbits an inner binary on a much wider orbit. A hierarchical 
triple system can be treated as two binaries on slowly perturbed Kepler orbits by using Jacobi 
coordinates, where the position of the secondary of mass mi of the inner binary is relative to the 
primary of mass mo and the position of the third body of mass m2 is relative to the center of mass 
of mo and mi. A hierarchical two-planet system is simply a hierarchical triple system with mi and 
m2 much smaller than mo- In § 2 we derive the orbital parameters in Jacobi coordinates obtained 
by the observers from a two (or more generally multiple) Kepler fit to the radial velocity variations 
of a host star and show that Jacobi coordinates should be used in any analyses of hierarchical (and 
possibly other types of) planetary systems. Star-centered or astrocentric coordinates can introduce 
significant high-frequency variations in orbital elements that should be nearly constant on orbital 
timescales. The high-frequency variations can then lead to erroneous sensitivity of the evolution of 
the orbital elements to the starting epoch. 

Secular perturbation theories based on an expansion in a = 01/02 have been developed for 
hierarchical triple systems. The expansion to order a^, called the quadrupole approximation, was 
developed by Kozai (1962) and Harrington (1968). This quadrupole-level secular perturbation the- 
ory is successful in explaining the Kozai mechanism, whereby the perturbations between the inner 
binary and the third body can lead to large variations in the eccentricity of the inner binary and 
the mutual inclination angle between the inner and outer binaries if the initial mutual inclination is 
sufficiently high. However, the quadrupole approximation is not adequate for studying hierarchical 
two-planet systems. Although there is no direct information on the mutual inclination angles in 
the known hierarchical two-planet systems, the formation of planets from a common disk of ma- 
terials surrounding the host star makes nearly coplanar orbits the most probable configuration. 
When the orbits are coplanar, the conservation of total angular momentum of the system and the 
secularly constant semimajor axes mean that the eccentricities of the two orbits are coupled and 
oscillate out of phase. The quadrupole term does not contribute to these eccentricity oscillations 
for coplanar orbits. Marchal (1990), Krymolowski & Mazeh (1999), and Ford, Kozinsky, & Rasio 
(2000) have extended the approximation to octupole (i.e., a^) order. Blaes, Lee, &: Socrates (2002) 
have recently applied this octupole-level secular perturbation theory, with a sign error corrected 
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(see footnote 5 below) and with modifications to include the effects of general relativistic precession 
and gravitational radiation on the inner binary, to study the dynamical evolution of hierarchical 
triples of supermassivc black holes. Unlike the quadrupole term, the octupole term does produce 
eccentricity oscillations for coplanar orbits. 

In this paper we use both the octupolc-level secular perturbation theory and direct numerical 
orbit integrations to investigate the dynamical evolution of coplanar hierarchical two-planet sys- 
tems. The applicability of the octupole theory is limited by its use of an expansion in a and an 
averaging over the inner and outer orbital motions. In particular, the orbit averaging eliminates 
the effects of mean-motion commensurabilities and the possible development of instabilities. We 
establish the validity and limits of the octupole theory by comparison with direct numerical orbit 
integrations. In § 3 we summarize the derivation of the octupole theory, compare the octupole 
theory to the classical Laplace-Lagrange secular solution, and deduce some useful results from the 
octupole equations analytically. In § 4 we study the dynamical evolution of the HD 168443 and 
HD 12661 systems and their variants. The HD 168443 system is not in a secular resonance and 
its secular resonance variable w\ — W2 circulates, where w\^2 are the longitudes of periapse of the 
inner and outer orbits, respectively. For the HD 168443 system and a wide variety of systems 
with a ~ 0.1 (including some for which the octupole theory predicts rather unusual dynamical 
behaviors), we show that the octupole results are in excellent agreement with the direct integration 
results. Direct integrations of two-planet systems similar to HD 168443, but with different initial 
02 (and hence different initial a), are used to show that systems with initial a above a critical 
value arc generally unstable. As anticipated by the analytic results derived in § 3, the HD 12661 
system is in a secular resonance with wi — W2 librating about 180°, and it shows large amplitude 
variations of both eccentricities. We show that the evolution of the HD 12661 system with the 
best-fit orbital parameters and sin i = 1 (i is the inclination of the orbital plane from the plane of 
the sky) is affected by the proximity to the 11:2 mean-motion commensur ability, but that small 
changes in the orbital period of the outer planet within the uncertainty can result in configurations 
that are not affected by mean-motion commensurabilities. For the latter type of configurations, we 
show that the octupole results are in reasonably good agreement with the direct integration results, 
even though a(w 0.32) is quite large. We also consider the effects of varying the inclination i and 
show that the HD 12661 system is unstable if sini < 0.3. Our conclusions are summarized in § 5. 



2. JACOBI ORBITAL PARAMETERS FROM MULTIPLE KEPLER FITS TO 

RADIAL VELOCITY OBSERVATIONS 

Except for systems such as GJ 876, where the perturbations between the two planets are 
significant on orbital timescales and a d5mamical fit to the stellar radial velocity variations is 

essential (Laughlin & Chambers 2001; Rivera & Lissauer 2001; Nauenberg 2002), it is often adequate 
to fit the radial velocity variations of a star with two (or more) planets over the time span of the 
available observations by assuming that the planets are on unperturbed Kepler orbits. Many 
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authors have assumed that the orbital parameters obtained from the multiple Kepler fits are in 
astrocentric coordinates, but Lissauer & Rivera (2001; sec also Note Added in Proof of Rivera 
& Lissauer 2000) have pointed out that Jacobi coordinates "better emulate" the assumption of 
unperturbed Kepler orbits. In this section we derive explicitly the orbital parameters in Jacobi 
coordinates obtained from a two-Kepler fit. It is straightforward to generalize the derivation to 
an A?^-Kepler fit. We show that especially for hierarchical systems such as HD 168443, the use of 
astrocentric coordinates can introduce erroneous features in the evolution of the orbital elements. 

Let us consider a system consisting of a central star of mass mo, an inner planet of mass mi, 
and an outer planet of mass mg, and use Jacobi coordinates, with ri being the position of mi 
relative to mo and r2 being the position of m2 relative to the center of mass of mo and mi. We 
shall refer to the orbit of mi relative to mo as the inner orbit and the orbit of m2 relative to the 
center of mass of mo and mi as the outer orbit. 

If the inner orbit is an unperturbed Kepler orbit, the line-of-sight (LOS) component of the 
velocity of mi relative to mo is 



l,r 



fi sin(a;i + /i) + ri/i cos(a;i + /i) 
— 27rai 



Pi0^ 

where a dot over a symbol denotes d/dt and ai, ei, ii, wi, /i, and 



smzi 

[cos(a;i + /i) + ei coswi] sinii, (1) 



(2) 



A/G(mo + mi) 

are, respectively, the semimajor axis, eccentricity, inclination, argument of periapse, true anomaly, 
and period of the inner orbit. In equation (2) G is the gravitational constant. Note that the 
reference plane is the plane of the sky, and the planet approaches the observer at the ascending 
node, but the radial velocity of an approaching planet is negative. Then the LOS component of 
the velocity of mo relative to the center of mass of mo and mi is 

Vo r = ; yi,r = Ki [cos(a;i + /i) + ei cos ui] , (3) 

mo + mi 



where the amplitude 

(mo + mi)2/3 0:3^- 



Similarly, if the outer orbit is an unperturbed Kepler orbit, the LOS component of the velocity 
of m2 relative to the center of mass of mo and mi is 

^2,r = „ o [C0s(a;2 + /2) + 62 COS UJ2] sin 12, (5) 

-P2VI -^2 
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where 02, 62, ^2, a;2, /2, and 



27ra' 



3/2 



P2 = " (6) 

^/G{mo + mi+ 

are, respectively, the semimajor axis, eccentricity, inclination, argument of periapse, true anomaly, 
and period of the outer orbit. Then the LOS component of the velocity of the center of mass of mo 
and mi relative to the center of mass of the whole system is 

^01 r = ^2,r = K2 [cOs{l02 + /2) + 62 COSa;2] , (7) 

mo + mi + m2 



where the amplitude 



K - ( ^^^ ^ "T-2 sini2 1 , . 

^~\P2) (mo + mi + m2)2/3 0^- ^> 



Thus, if the orbits of the planets in Jacobi coordinates are unperturbed Kepler orbits, the 
radial velocity of the star mo is 

Vr = ^oV + ^01,.+^ 

= Ki [cos(a;i + /i) + ei cos ui] + K2 [cos(a;2 + /2) + 62 cos 0J2] + V, (9) 

where V is the LOS velocity of the center of mass of the whole system relative to the observer. 
Equation (9) is exactly the formula used by observers in two-Kepler fits, but with the amplitudes 
defined in equations (4) and (8) and the orbital periods defined in equations (2) and (6).^ Since the 
true anomaly fj depends on Pj, ej, and the time of periapse passage Tp^^j-ij, a two-Kepler fit directly 
yields five parameters Pj, Kj, Cj, cvj, and Tpcrij for each orbit. Table 1 shows these parameters 
for the HD 168443 planets from Marcy et al. (2001) and the HD 12661 planets from D. A. Fischer 
(2002, private communication). If the stellar mass mo is known, equations (2), (4), (6), and (8) can 
be used to derive rrij and aj for assumed sin ij . (A common but less accurate practice is to derive 
rrij sinzj and aj by assuming that mi and m2 are negligible compared to mo.) Table 1 also shows 
ruj (in units of Jupiter mass Mj) and aj for sini^ = 1. Throughout this paper, we refer to the 
substcUar companions of HD 168443 as planets, but it should be noted that their sin = 1 (i.e., 
minimum) masses are about 7.7 and 17Mj, which are near or above the deuterium-burning limit. 

The alternative assumption that the orbits of the planets in astrocentric coordinates are unper- 
turbed Kepler orbits would also yield equation (9) for the radial velocity of the host star, but with 
Kj = [2TTG{mo+mj)/Pj]^/^mjsmij{mo + mi+m2y^il-e'^j)-'^/^ and = 2Traf'^[Gimo+mj)]-^/'^ . 
However, as we show next, while the Jacobi orbits of the planets of a hierarchical system are nearly 
Keplerian on orbital timescales, the astrocentric orbit of the outer planet can deviate significantly 
from a Kepler orbit on orbital timescales. 



^Note that the planetary masses in the equations for the amplitudes and the orbital periods are the physical 
masses rrij and not the Jacobi masses mj Xife=o '^fe/ Xife=o ^ stated by Lissauer & Rivera (2001). 
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Table 1. Orbital Parameters of the HD 168443 and HD 12661 Planets 



Parameter 


HD 168443 


HD 12661 


Inner 


Outer 


Inner 


Outer 


P (days) 


58.10 


1770 


263.3 


1444.5 


K (ms-^) 


472.7 


289 


74.4 


27.4 


e 


0.53 


0.20 


0.35 


0.20 


(J (deg) 


172.9 


62.9 


292.6 


147.0 


Tperi (JD) 


2450047.58 


2450250.6 


2449943.7 


2449673.9 


m [Mj) 


7.73 


17.23 


2.30 


1.57 


a (AU) 


0.295 


2.90 


0.823 


2.56 



Note. — The parameters P, K, e, lo, and Tpdi from two-Kepler 
fits by Marcy et al. 2001 and D. A. Fischer 2002, private commu- 
nication. The parameters m and a are derived for sini = 1 and 
adopted stellar masses of IMMq and 1.O7M0 for HD 168443 and 
HD 12661, respectively. 
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(a) HD 168443, sinj = l, astrocentric coords 
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(b) HD 168443, sini = l, Jacobi coords 
I 1 I ^ I , I 1 I , I , 



H 1 1 \ 1 \ 1 1 1 \ H 



Octupole 

H 1 1 1 1 1 1 1 1 1 1 1 1 h- 




5x10* 10" 
t (yr) 



1.5x10" 



Fig. 1. — (a) Variations in the orbital semimajor axes, ai and a2, and eccentricities, ei and 62, of 
the HD 168443 planets with sini = 1, if we assume that the best-fit orbital parameters obtained 
from the two-Kepler fit are in astrocentric coordinates and plot the astrocentric aj and ej. The 
solid and dotted lines are from direct numerical orbit integrations starting at Tperi,i and rperi,2 of 
Table 1, respectively. (6) Same as (a), but we interpret the best-fit orbital parameters obtained 
from the two-Kepler fit as orbital parameters in Jacobi coordinates and plot the Jacobi aj and ej. 
The solid and dotted lines are almost indistinguishable. The dashed lines in the lower two panels 
of (6) are from the octupole-level secular perturbation theory. 



Figure la shows the variations in the semimajor axes and eccentricities of the HD 168443 
planets if we assume that the orbital parameters obtained from the two-Kepler fit are in astrocentric 
coordinates and plot the astrocentric aj and ej. The planets are assumed to be on coplanar 
orbits with sinz = 1 (note that coplanar orbits have the same inclinations, ii = 12 = i, and the 
same longitudes of ascending node, = ^1, referenced to the plane of the sky). The solid and 
dotted lines are from direct numerical orbit integrations starting at Tj^cj.; 1 and TpcTi.2 

of Table 1, 

respectively. The direct integrations were performed using the Wisdom-Holman (1991) integrator 
contained in the SWIFT^ software package, with input and output in astrocentric coordinates. 
Figure 16 is similar to Figure la, but it shows the variations in the Jacobi aj and ej, with the 
initial orbital parameters also in Jacobi coordinates. The direct integrations were performed using 
the modified Wisdom-Holman integrator described in § 4. In astrocentric coordinates (Fig. 1 a) , the 



''See http://www.boulder.swri.edu/~lial/swift.litml. 
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evolution of the orbital elements is sensitive to the starting epoch; and while ai is nearly constant 
and ei oscillates only on the secular timescalc, both a2 and 62 show significant fluctuations on 
the orbital timescales (note that we have reduced overcrowding in Fig. la by using a sampling 
interval of 50 yr, which is long compared to the orbital timescales). Plots of the variations in the 
orbital elements of the HD 168443 planets in Marcy et al. (2001), Nagasawa, Lin, k. Ida (2003), 
and Udry, Mayor, k. Queloz (2003, from calculations by W. Benz) show similar behaviors and are 
probably due to the use of astroccntric coordinates. In contrast, in Jacobi coordinates (Fig. 16), 
the evolution is not sensitive to the starting epoch (the solid and dotted lines in Fig. lb are almost 
indistinguishable); both oi and 02 are nearly constant, while both ei and 62 oscillate only on the 
secular timescale, with the maximum in ei coinciding with the minimum in 62 and vice versa. As 
we shall see, these behaviors can be understood with the octupole-level secular perturbation theory 
for coplanar hierarchical two-planet systems, which is also based on Jacobi coordinates. 

Prom the facts that the Jacobi aj are nearly constant and that the Jacobi Cj oscillate only on 
the secular timescale for hierarchical systems, it is easy to estimate the fractional fluctuation in 
the astrocentric 02, which is primarily due to the velocity of the star mo relative to the center of 
mass of mo and mi. It is approximately 4(mi/mo)[(a2/fi'i)(l + ei,max)/(l — ei^max)]^^^ when the 
Jacobi ei is at its maximum ei^max (and 62 is at its minimum, which is assumed to be small and 
is neglected). The fluctuations in the astrocentric 02 in Figure la are in good agreement with this 
estimate. Note that the fractional fluctuation is larger for smaller 01/02, i.e., for systems that are 
more hierarchical. 

In their study of the three-planet v And system, Rivera k, Lissauer (2000) have also reported 
that the use of Jacobi coordinates eliminates the high-frequency variations in the semimajor axes 
and eccentricities of the outer two planets and reduces the sensitivity of the evolution to the initial 
epoch. It is clear that Jacobi coordinates should be used in the analysis of hierarchical systems 
such as HD 168443 (where 01/02 ~ 0.10) and v And (where 01/02 ~ 0.071 and 02/03 ~ 0.33). 
It is likely that Jacobi coordinates should also be used in the analysis of other types of planetary 
systems, especially those for which multiple Kepler fits are adequate. 

3. OCTUPOLE-LEVEL SECULAR PERTURBATION THEORY 

As we mentioned in § 1, an approximate theory that describes the secular evolution of copla- 
nar hierarchical two-planet systems in Jacobi coordinates, such as that shown in Figure Ih for the 
HD 168443 system, is the octupole-level secular perturbation theory. In this section we summa- 
rize the derivation of this octupole theory, compare it to the classical Laplace-Lagrange secular 
perturbation theory, and deduce some results from the octupole equations analytically. 
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3.1. Equations 

As in §2, we consider a system consisting of a central star of mass mo, an inner planet of mass 
mi, and an outer planet of mass m2, and use Jacobi coordinates, with ri being the position of mi 
relative to mo and r2 being the position of m2 relative to the center of mass of mo and mi. The 
Hamiltonian of this system is 

^ ^ _Gmomi _ G{n->o + nn)m2 _ ^^^^^ ( ± - ^) - Gmim, (---], (10) 



2ai 2a2 \ro2 r2 J \ri2 ri 

where is the osculating semimajor axis of the jth orbit (with j = 1 and 2 for the inner and outer 
orbits, respectively), and r\a is the distance between m^ and m2. With r\ < ri, both l/ro2 and 
l/ri2 can be expanded in powers of ri/ri, leading to the expanded Hamiltonian 



GmQmi G(mo + mi)m2 G 

rC = — 



--^^a^Mj^) (^) P,(cos$), (11) 

where a = ai/ai. 



M, = momim2^^^V-^T^' (^2) 
(mo + mi)'= 

Pfe is the Legendre polynomial of degree k, and $ is the angle between ri and ri- The first 
two terms in equations (10) and (11) represent the independent Kepler motions of the inner and 
outer orbits, while the remaining terms represent the perturbations to the Kepler motions. For 
hierarchical systems with ri <C r2, the Jacobi decomposition leads to two slowly perturbed Kepler 
orbits, even if mi and m2 are not much smaller than mo. 

In the general case where mutually inclined orbits are allowed, it is convenient to use the 
Delaunay variables with the invariable plane as the reference plane. The coordinates are the mean 
anomalies Ij, the arguments of periapse gj = LOj, and the longitudes of ascending node hj = ^j, 
and the conjugate momenta are 



Li = "^""^^ V G{mo + mi)rai , (13) 
mo + mi 

(mo + mi)m2 rr^, ■ ■ r — 

L2 = \/G{mo + mi + m2)a2, (14) 

mo + mi + m2 

Gj = L.^Jl^^, (15) 
Hj = Gj cos ij, (16) 

where ej and ij are the eccentricities and inclinations of the orbits. The momenta Lj, Gj, and Hj 
are, respectively, the magnitude of the maximum possible angular momentum (if the orbit were 
circular), the magnitude of the angular momentum, and the z-component of the angular momentum 
of the jth orbit. By expressing cos $ = ri-r2/(rir2) in terms of the Delaunay variables, it is easy to 
show that the Hamiltonian contains hi and h2 only in the combination hi — hi, and hi — h2 = 180° 
if the invariable plane is used as the reference plane. Constant hi — hi and constant total angular 
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momentum Hi + H2 mean hj and Hj can be eliminated from the problem — usually referred to 
as the elimination of nodes — and the system is reduced from 6 to 4 degrees of freedom (see, e.g., 
Marchal 1990). 

Kozai (1962) studied the secular evolution of hierarchical triple systems with mi <^ m-2 <C rriQ 
and 62 = by retaining just the lowest-order (a^) quadrupole term in the series in equation (11) 
and averaging the Hamiltonian over the inner and outer orbital motions. Harrington (1968) gener- 
alized the quadrupole analysis to general hierarchical triple systems. Marchal (1990), Krymolowski 
&; Mazeh (1999), and Ford et al. (2000) extended the approximation by retaining also the oc- 
tupole term of order in the series in equation (11) and derived an octupole-level orbit-averaged 
Hamiltonian. The approach used in these studies to average the Hamiltonian is the von Zeipel 
method, which involves the determination of a canonical transformation such that the transformed 
Hamiltonian is independent of the transformed mean anomalies li and I2, and such that the original 
and transformed variables differ only by high frequency terms that are multiplied by powers of the 
small parameter a. Since the transformed Hamiltonian is independent of the transformed Zj, the 
transformed Lj and the corresponding semimajor axes are constant. 

In this paper we focus on hierarchical two-planet systems with coplanar orbits, as the formation 
of planets from a common disk of materials surrounding the host star makes nearly coplanar orbits 
the most probable configuration. In the limit of coplanar orbits {i\ = 12 = 0°), the Delaunay 
variables described above are not well-defined, since there are no ascending nodes with respect 
to the invariable plane, but the longitudes of periapse Wj are well-defined. Thus, for coplanar 
orbits, it is convenient to use the canonical variables lj, Wj, Lj, and Gj. In the limit of coplanar 
orbits, by noting that the expression for the angle if between the directions of periapse reduces to 
cosip = — cos(5i — 52) = cos{gi — 52 + 180°) = cos(iui — ^2) (since h\ — h2 = 180°), the doubly 
averaged octupole-level Hamiltonian derived by Marchal (1990), Krymolowski &; Mazeh (1999), 
and Ford et al. (2000) can be written as^ 

^oct = -77} ^° \t2 - o( I 2 ,;2 -2C2 2 + 36^ +^36162 4+3ef COS tUi -tI72 , 17 

2(mo+mi)L| 2(mo + mi + 7^2)^2 



where 



C2 = 1 G'^jmo + rmYml Lj 

16 (mo + mi + m2)^(momi)^ L2G2 ' 

_ 15 G'^(mo + mi)'^m|(mo - mi) Lf 
^ ~ 64 (mo + mi + m2)^(momi)5 L^G^' ^ ' 



■''Thoro is a sign error common to both Krymolowski & Mazeh (1999) and Ford et aL (2000). We follow Ford et al. 
and denote the coefficients of the quadrupole and octupole terms by C2 and C3, respectively. Krymolowski & Mazeh 
denote these coefficients by Ci and C2, respectively. In the averaged Hamiltonian and the subsequent equations of 
motion of Ford ct al., the coefficient C3 should be replaced by — C3. Similarly, in the averaged Hamiltonian and the 
subsequent equations of motion of Krymolowski & Mazeh, the coefficient C2 should be replaced by — C2. 
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and 

(20) 

As in Ford et al. (2000), we do not include in equation (17) terms of order a^/^ induced by the 
canonical transformation of the von Zeipel method, which were included partially by Marchal (1990) 
and fully by Krymolowski & Mazeh (1999). 

As we mentioned above, the orbit averaging eliminates Ij from TYoct) and hence Lj (and aj) are 
constant. Furthermore, Hoct contains wi and W2 only in the combination wi — W2, and choosing 
this difference as a new variable in a canonical transformation simply reveals that the total angular 
momentum Gi + G2 is an integral of motion. Thus, in the octupole approximation, the coplanar 
problem is reduced to one degree of freedom with the last two terms in equation (17) as an integral 
of motion. We can also see from equation (17) that the quadrupole approximation is not adequate 
for studying coplanar systems. If we drop the octupole term with coefficient C3 in equation (17), the 
Hamiltonian is independent of wi and W2, meaning that both ei and 62 are constant, which is not 
generally true. (For mutually inclined orbits, only 62 is constant in the quadrupole approximation; 
Harrington 1968.) 

From the equations of motion 

dujj _ dHoct , dOj _ dUoct 

dt ~ dC, ^"""^ dt ~ dwj ' 

we obtain the following equations for the variation of ej and zuj: 
de, . (1 - )V^ (1 + f ef ) . 

-df = -^^^^^ (TT^^ sm(..,-..2), (22) 

^ = A2,e,^^-^-^sm{w,-W2), (23) 
dt (1 - eiY 

d^i , (1 - ef )V^ , fe2\ {l-eir^ (l + fef) 

- ^117; - ^12 — 7--, 2,^/0 cos(roi - W2), (24) 



dt ^^(i_e2)3/2 ''\ei) (l-e2)5/2 

d^2 , (l + jef) , /eA (l + 4e|)(l+|e|) 

- ^22 Ti Srf - -421 — -r-^ — cos(tx7i - W2), (25) 



dt ''\e2j (l-ei)3 

where the constant coefficients 



4:C3 f G2\^ 15 f m2 \ f mo - mi\ ^ 



Li \L2 J 16 \mo + mi/ \mo + mi/ 

4C3/G2y 15 momi (niQ-mi^ 3 

Ml = — — = —n2- ^ a\ 26b 

L2 \L2 J 16 (rriQ + rniY \mQ + mi' 
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A 



11 



A 



22 



12C2 

VIC2 
L2 



G2 
L2 

G2 
L2 



3 

-n2 



m2 



mo + mi 
moiTii 



a 



4 (mo + mi)^ 



a 



with a = 01/02 and the mean motions 

m = [G(mo + mi)/a5]^^^ 



and n2 = [G(mo + mi + m2)/a2] 



31 1/2 



(26c) 
(26d) 

(27) 



Consistent with the fact that the system can be reduced to one degree of freedom, either equation 
(22) or (23) is redundant, since ei and 62 arc related by the conservation of total angular momentum, 
and equations (24) and (25) can be combined as 



d{wi - W2) 
di 



A 



11 



(1 - ef )V2 fL,^ (1 + |ef ) 
L2, 



(1 



ei)3/2 



(1 



=2N2 



112 



e2 ^ (1 - e?)V2 (1 + |e2) 
ei. 



(1 - ei)V2 

X COs(tI7l — tI72). 



L2 



ei\ (l + 4ei)(l + fef) 
62 ; (l-ei)3 



(28) 



If we relax the assumption of coplanar orbits and allow a small mutual inclination angle imu 
between the inner and outer orbits, the number of degrees of freedom is increased from one to two, 
but an expansion of the orbit-averaged octupole-level Hamiltonian of Marchal (1990), Krymolowski 
& Mazeh (1999), and Ford et al. (2000) in powers of imu shows that equations (22)-(25) are only 
modified by terms of order and higher. One of the additional terms for dei/dt (which has 
an octupole term only in the coplanar limit; see eq. [22]) is a quadrupole term with coefficient 
C2^mu- Thus we expect a small mutual inclination angle Zmu ^ "^^^ to have little effect on the 
secular evolution of ej and zuj of a hierarchical system. This is confirmed by direct numerical orbit 
integrations in § 4. 

It is important to note that the octupole-level secular perturbation theory is derived under 
the assumptions that there are no mean-motion commensurabilities and that a (or more precisely 
^1/^2) ^1- As we shall sec in § 4, close proximity to even a rather high order commensurability 
(e.g., 11:2) can affect the evolution of a system. We shall also see in § 4 that although there is a 
limit to how large a can be for a system to be stable, the octupole theory can provide a reasonable 
description of the secular evolution for a as large as 1/3 if the planets are not too massive. 



3.2. Comparison with Classical Laplace-Lagrange Secular Perturbation Theory 

As we mentioned in § 1, the classical secular perturbation theory developed for the solar system 
is based on an expansion in the eccentricities and inclinations. While it is valid to all orders in 
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the ratio of the semimajor axes a, it can be applied only to systems with small planetary masses 
on nearly circular and nearly coplanar orbits. In contrast, the octupole-levcl secular perturbation 
theory is based on an expansion in a, and it can be applied to hierarchical systems with a wide 
range of masses, eccentricities and, in its most general form, inclinations. We now show that the 
octupole equations (22)-(25) agree with the equations for the classical Laplace-Lagrange secular 
solution for coplanar systems where the planetary masses, eccentricities, and a are all small. 



To the lowest order in the eccentricities, equations (22)-(25) reduce to 

-^1262 sin(roi — W2), (29) 



dei 

'dt 
de2 

It 

An - Aui — ] cos(roi - ZU2), (31) 



^2iei sin(roi - UJ2), (30) 



at \ei 

"^"^^ A22 - A21 (- ) cos(ti7i - W2), (32) 



dt \e2 

where A12 = (15/16)(m2/mo)nia'^, A21 = (15/16)(mi/mo)n2Q!^, An = (3/4)(m2/mo)nia^, and 
^22 = (3/4)(mi/mo)n2a^ in the limit mi,m2 <C mo. 

The classical Laplace-Lagrange secular solution is based on retaining just the secular terms in 
the disturbing functions up to second order in the eccentricities and inclinations. For a coplanar 
planetary system, the equations for the variation of ej and Wj (see, e.g., Murray &; Dermott 1999) 
are of the same form as equations (29)-(32), but with the Ajk replaced by 

A' "^2 2,(2) /X /OQ N 

^2 = 4"^ mo + mi '' ^3/2(«)' (33a) 
= \n2 ab^^Ua), (33b) 

A[i = ]ni ""^ a¥'ha), (33c) 

^22 = T«2 — ^ — a4/2(«)' (33d) 
4 mo + m2 ' 

where b^J^{a) is the Laplace coefficient. (It should be noted that the classical theory uses astro- 
centric coordinates instead of Jacobi coordinates, but that the distinction between the two sets 
of coordinates vanishes in the limit mi/rriQ <^ a^/^/4; see §2.) Since h^^j^{a) = (15/4)a^ and 

^3/2 ("^-^ = 3a to the lowest order in a, the classical and octupole secular perturbation equations are 
identical in the limit of small planetary masses, eccentricities, and a. 
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3.3. Some Analytic Results 



Although it does not appear that the octupole-level secular perturbation equations can be 
solved analytically, some useful results can be deduced analytically. We begin by rewriting equations 
(22), (23), and (28) as 



de2 



= -Pe2 



pXei 



(1 - ef )V^ (1 + f ef ) 
(1 - ei)V2 

(1 + hi) 



sin(roi —W2), 



elf 



sin(cc7i — -072), 



(1 

(l-ef)V^ 

(1 - ei)3/2 " (1 

" 62 \ (1 - e?)V2 (1 + fef) 



(34) 



(35) 



:i + i^f ) 

=2\2 



ei 



X cos(cc7i 



(1 -e: 



2)5/2 



-A 



(l + 4ei) (l + |ef) 
(1 - elf 



where 



A 

r 



and 



A12 


A21 


All 


A22 


A22 


A21 


All 


A12 






1 


4 



5 / mo — mi 



4 ymo + mi 
L2' 



a. 



11 



3a3 



mo + mi 
m2 



1 



(36) 

(37) 

(38) 
(39) 

(40) 



Recall that either e\ or 62 can be eliminated from the above equations by the conservation of total 
angular momentum, which can be written in dimensionless form as 7 = {Gi + G2)/{Li + L2) = 
[A(l - ef )V2 + (1 _ ei)V2]/(A + 1) = constant. 

We can see from equations (34)-(36) that coplanar hierarchical systems with the same /3 and 
A (which are constant) and the same initial ei, 62, and zui — tU2 have the same trajectory in the 
phase-space diagram of ei (or 62) versus wi — W2 and can differ only in the period of eccentricity 
oscillations. In particular, in the limit mi, 1712 <C mo, since /3 ~ (5/4)a, A ~ (mi/m2)Q!^/^, and 
te ~ 4(mo/m2)/(3a3ni), systems with the same a, mi/m2, and initial ei, 62 and wi — W2, but 
different mo, m2 and/or ai, differ only in the period of eccentricity oscillations, which is proportional 
to (mo/m2)/ni. As we discussed in § 2, a two-Kepler fit yields Cj, Uj, and in the limit mi, m2 ^ mo, 
rrij sinij and Oj, but not the inclinations ij of the orbital planes to the plane of the sky. (Hereafter, 
all inclinations are relative to the plane of the sky and not to the invariable plane.) If we assume 
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that the orbits are coplanar, the above results mean that the trajectory in ei (or 62) versus wi — W2 
and the amphtudes of eccentricity oscillations should be independent of the unknown inclination i 
of both orbits, and that the period of eccentricity oscillations should be proportional to sini. 

The phase-space structure of the secular evolution of coplanar hierarchical two-planet systems 
can be understood by plotting the trajectories of systems with the same /?, A, and 7 in a diagram 
of ei (or 62) versus wi — W2 (sec, e.g., Figs. 2 and 6 below). The fixed points in the phase-space 
diagram can be found by solving dei/dr = de2/dT = d{wi — W2)/dT = 0. It is clear from equations 
(34) and (35) that dei/dr = de2/dT = requires sin(tA7i — W2) = or oji — ^2 = 0° or 180°. As 
we shall see, there are usually two libration islands, one about a fixed point at wi — ZU2 = 0° and 
another about a fixed point at wi — W2 = 180°; additional fixed points are possible if the total 
angular momentum is low. 

At wi — 'CU2 ~ 90° and 270°, since cos(ti7i — ZU2) ~ 0, equation (36) reduces to d{wi—W2)/dT ~ 
(1 — ef )^/^/(l — 62)^/^ — A(l -I- 3ef/2)/(l — 62)^. Therefore, to the lowest order in the eccentricities, 
if A 1, d{zui — W2)/dT « and wi — W2 should be nearly constant while both ei and 62 change. 
According to equations (34) and (35), ei should be decreasing (from near its maximum possible 
value for the given total angular momentum to near zero) and 62 should be increasing (from near 
zero to near its maximum possible value for the given total angular momentum) at vj\—W2^ 90° , 
where sin(cc7i — W2) ~ 1, and vice versa at w\ — W2 ~ 270°, where sin(ci7i — ■072) ~ —1. These 
behaviors are most consistent with a phase space that is dominated by two large libration islands, 
one about a fixed point at vj\ — W2 = 0° and another about a fixed point at wi — W2 = 180°. 
Note also that for a given total angular momentum the maximum possible values of ei and 62 
are comparable if A pa 1. Finite eccentricities change the condition for d{wi —W2)/dT k, at 
■uji-vj2^ 90° and 270° to A ?a (1 - ef )V2(i _ elf/'^/il + 3ef/2), which is slightly less than 1 
for moderate eccentricities. To estimate how much smaller than 1 the critical value of A is for a 
given dimcnsionless total angular momentum 7, wc substitute into the above condition the value 
of the eccentricities when they are equal [ei = 62 = (1 — T'^)^^^] and obtain Acrit = ^7^/(5 — 37^). 
Note that Acrit = 1 is recovered for 7 = 1. Therefore, large libration islands and large amplitude 
variations of both ej are likely if A Acrit = 27^/(5 — 37^). 

4. NUMERICAL RESULTS 

In this section we study the dynamical evolution of the HD 168443 and HD 12661 systems and 
their variants and demonstrate the validity and limits of the octupole-level secular perturbation 
theory by comparison with direct numerical orbit integrations. Except for the direct integrations 
discussed in § 4.3, the planets are assumed to be on coplanar orbits. Since the octupole-level 
secular perturbation equations do not involve high frequency terms, numerical integrations of these 
equations are rapid. We integrated equations (35) and (36), with ei found from conservation of 
total angular momentum, using a Bulirsch-Stoer integrator. The direct numerical orbit integrations 
were performed using a modified version of the Wisdom-Holman (1991) integrator contained in the 
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SWIFT software package. In addition to changing the input and output to Jacobi orbital elements, 
we divide the Hamiltonian into a part that describes the Kepler motions of the inner and outer 
orbits and a part that describes the perturbations to the Kepler motions using equation (10) instead 
of the division used by Wisdom & Holman (1991), which moves the term Gmim2/r2 from the 
perturbation Hamiltonian to the Kepler Hamiltonian of the outer orbit. This modified integrator 
can handle hierarchical systems where mi and m2 are not small. 

4.1. HD 168443 

The best-fit orbital parameters of the HD 168443 planets from Marcy et al. (2001) and the 
inferred planetary masses and semimajor axes for sini = 1 are listed in Tabic 1. This system 
has a = ai/a2 = 0.102, /3 = 5(mo - mi)a/[4(mo + mi)] = 0.126, A = L1/L2 = 0.143, and 
7 = (Gi + G2)/{Li + L2) = 0.963. Figure 16 shows the variations in the Jacobi aj and Cj, with the 
solid and dotted lines (almost indistinguishable) from direct numerical orbit integrations starting 
at Tperi,i and rperi,2 of Table 1, respectively, and the dashed lines in the lower two panels from 
the octupole-level secular perturbation theory. The trajectories in Figure 2 marked by the squares 
(which denote the current parameters of the HD 168443 system) are the trajectories of the sinz = 1 
HD 168443 system in the phase-space diagrams of ei versus wi — W2 (upper panels) and 62 versus 
wi — VJ2 (lower panels), with the left and right panels showing the results from direct integrations 
and the octupole theory, respectively (see below for a discussion of the other trajectories in Fig. 2). 
The direct integration results agree with the octupole theory in that the evolution is not sensitive 
to the starting epoch and that ai and 02 are nearly constant. In addition, the trajectories in the 
diagram of ei (or 62) versus wi — W2 from the direct integrations are in excellent agreement with 
that from the octupole theory. The only noticeable difference is that the period of eccentricity 
oscillations predicted by the octupole theory is about 3% longer that that found in the direct 
integrations (« 1.8 x 10^ yr). Krymolowski & Mazeh (1999) have reported that the oscillation 
period predicted by the octupole theory can be improved by the inclusion of the terms of order 
a^l"^ induced by the canonical transformation of the von Zeipel method and neglected by us. 

To study the effects of sinz, we have performed two direct integrations of the coplanar HD 
168443 system with sinz = 0.4, one starting at Tperi,! and the other rperi,2- (Although the two- 
Kepler fit to the radial velocity observations does not yield sinij, the lack of evidence for stellar 
wobble in the Hipparcos astrometric data limits sini2 ^ 0.4 for the outer planet of HD 168443; 
Marcy et al. 2001.) In both cases, the amplitudes of eccentricity oscillations and the trajectory 
in ei (or 62) versus tui — W2 are almost identical to those shown in Figures lb and 2 for sini = 1, 
and the factor by which the eccentricity oscillation period shortens agrees with sini = 0.4 to better 
than ±0.001. These results are in good agreement with the analytic results derived in § 3.3 in the 
limit mi,m2 <S mo, even though m2 is almost 4AM j (and 1712/ tuq 0.042) when smi = 0.4. 

In addition to the trajectories of the sini = 1 HD 168443 system. Figure 2 also shows the 
trajectories of systems with the same masses, initial semimajor axes, and total angular momentum 




Fig. 2. — Trajectories in the phase-space diagrams of ei vs. wi — W2 {upper panels) and 62 vs. 
wi — VJ2 {lower panels) for two-planet systems with the same masses, initial semimajor axes, and 
total angular momentum as the sini = 1 HD 168443 system. The left and right panels show 
the results from direct numerical orbit integrations and the octupole-level secular perturbation 
theory, respectively. The trajectories through the squares are those of the HD 168443 system, 
with the squares showing the current parameters of the system. The initial conditions for the 
other trajectories are tui - ^2 = 0° and d = 0.06, 0.08, 0.10, 0.15, 0.25, 0.40, 0.65, or 0.70, 
or wi — VJ2 = 180° and e\ = 0.706 or 0.7072, with 62 being determined from the total angular 
momentum. All the direct integrations (except that for the HD 168443 system) start with the 
outer planet at apoapse and the inner planet at opposition. 



as the sini = 1 HD 168443 system. (Recall that the octupole results are also valid for other 
systems with the same (5, A, and 7.) There is excellent agreement between the direct-integration 
and octupole results in all cases. The HD 168443 system is not in a secular resonance and its 
w\ — W2 circulates. Indeed it is far from the relatively small libration islands about the fixed 

points at {vji —TU2,ei) = (0°, 0.046) and (180°, 0.702) in the phase-space diagram of ei versus 
wi — W2- The small libration islands and modest eccentricity variations in Figure 2 are consistent 
with the fact that these systems have A = 0.143 far from Acrit = 0.836 for the dimensionless total 
angular momentum 7 = 0.963 of these systems. 

The above ei values for the fixed points were obtained from the octupole theory as the roots 
of ci(roi — W2)/dT = at roi — W2 = 0° and 180°, with d{wi — W2)/dT from equation (36) and 62 
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Fig. 3. — Values of ei for the fixed points at wi — W2 = 0° {left panels) and 180° {right panel) in 
the phase-space diagram of e\ vs. w\ — W2 (like Fig. 2) as a function of the dimensionless total 
angular momentum 7 = (Gi + G2)/{Li + L2), as determined by the octupole theory for systems 
with the same (3{= 0.126) and A(= 0.143) as the sini = 1 HD 168443 system. The dotted lines 
indicate the dimensionless total angular momentum 7 = 0.963 of the sini = 1 HD 168443 system, 
and the dashed lines indicate the maximum possible ei as a function of 7. See text for the meaning 
of the square in the upper left panel. 



from the conservation of total angular momentum. As we discussed in § 3.3, the octupole theory 
predicts that the fixed points must be at w\ — vji = 0° or 180° . To demonstrate how the octupole 
theory can be used to explore the parameter space rapidly, we have used the same root finding 
procedure to determine the number and positions of fixed points as a function of the dimensionless 
total angular momentum 7 = {Gi + G2)/{Li + L2) for systems with the same (3 and A as the 
sin i = 1 HD 168443 system. The e\ values for the fixed points at w\ — W2 = 0° and 180° as a 
function of 7 are shown in the left and right panels of Figure 3, respectively. For all values of 7 
shown in Figure 3, there is an elliptic fixed point at wi — vj2 = 0° with relatively small e\ (lower 
left panel) and an elliptic fixed point at wi — ^2 = 180° with ei close to the maximum possible 
value indicated by the dashed line (right panel). These two fixed points are the ones seen above 
for systems with the same 7(= 0.963; dotted lines in Fig. 3) as the sin z = 1 HD 168443 system. 

For 0.872 < 7 < 0.8818, there are two additional fixed points at wi — W2 = 0° with ei 
very close to 1 (upper left panel of Fig. 3). These additional fixed points emerge at 7 0.8818 
with the same ei value (square in the upper left panel of Fig. 3), and the one with ei value that 
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increases (decreases) with decreasing 7 is an elliptic (hyperbolic) fixed point. Direct numerical 
orbit integrations confirm that there are indeed stable librations about an elliptic fixed point close 
to the additional one predicted by the octupole theory. For example, for 7 = 0.88, the octupole 
theory predicts that the additional elliptic fixed point is at ei = 0.9948 (and 62 = 0.1302). A direct 
integration of a system with these values of ei , 62 , and wi — W2 as initial conditions (and with the 
masses and semimajor axes of the sinz = 1 HD 168443 system and mean anomalies h = 0° and 
I2 = 180°) shows that this system has wi — wi librating about 0° with an amplitude of about 6° 
and d varying between 0.9946 and 0.9969. It should be noted that the treatment of this system 
with e\ > 0.9946 and ai ~ 0.3 AU as point masses interacting via Newtonian gravity is inadequate 
since the periapse distance of the inner planet is in fact less than the stellar radius. However, 
based on the octupole theory, we expect variants of this system with much larger ai (so that the 
periapse distance of the inner planet is well outside the stellar radius) but the same a to show 
similar libration behaviors. On the other hand, the inner planet could still come sufficiently close 
to the star for general relativistic precession and tidal effects to be important. For 7 < 0.872, the 
additional elliptic fixed point vanishes, leaving only the hyperbolic fixed point (upper left panel of 
Fig. 3). Both octupole and direct-integration calculations show that systems with 7 < 0.872, initial 
wi — W2 = 0° and initial ei above the hyperbolic fixed point have ei increasing to unity in a finite 
time. Again, we expect variants of these systems with larger ai (so that the periapse distance of 
the inner planet is well outside the stellar radius) but the same a to show similar increase in ei. 
However, before ei reaches unity, either the effects of tides and general relativistic precession would 
eventually stop the increase in ei, or the inner planet would collide with the star. 



4.2. HD 12661 

The best-fit orbital parameters of the HD 12661 planets from D. A. Fischer (2002, private 
communication) and the inferred planetary masses and semimajor axes for sinz = 1 are listed 
in Table 1. Figure 4 shows the variations in the Jacobi aj and ej, with the solid and dotted 
lines from direct numerical orbit integrations starting at Tperi,! and rperi,2 of Table 1, respectively, 
and the dashed lines in the lower two panels from the octupole-level secular perturbation theory. 
The direct integration results are not consistent with the secular theory in that the evolution of 
the orbital elements is sensitive to the starting epoch. It turns out that a system with exactly 
the best-fit orbital parameters is very close to the 11:2 mean-motion commensurability {P2/P1 = 
5.486 = 0.9975 x 11/2). The direct integration starting at Tperi,! (solid lines in Fig. 4) shows 
irregular fluctuations in aj and ej (note, e.g., the irregular jumps in the mean values of ai and 02 
at successive minima of 62) and is most likely chaotic. The direct integration starting at rperi,2 
(dotted lines in Fig. 4) does not show any obvious irregular jumps in aj or ej, and may be either 
regular or very weakly chaotic. But its smaller amplitudes of eccentricity variations are due to the 
11:2 commensurability. The chaos in one (and possibly both) of these calculations is due to the 
overlap of the resonances at the 11:2 commensurability (see Holman Sz Murray 1996 and Murray 
&; Holman 1997 for a similar situation in the planar elliptic restricted three-body problem). For 
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Fig. 4. — Variations in the Jacobi semimajor axes, ai and 02, and eccentricities, ei and 62, of 
the sini = 1 HD 12661 system with the best-fit orbital parameters hsted in Table 1. The solid 
and dotted lines are from direct integrations starting at Tperi,i and Tperi,2 of Table 1, respectively, 
and the dashed lines in the lower two panels are from the octupole theory. This system with 
initial P2/P1 = 0.9975 x 11/2 is affected by the 11:2 mean-motion commensurability, and the direct 
integration results are sensitive to the starting epoch. The direct integration starting at Tpgri,! 
shows irregular fluctuations in aj and ej (note, e.g., the irregular jumps in the mean values of a\ 
and 02 at successive minima of 62) and is most likely chaotic. The most noticeable effect of the 11:2 
commensurability on the direct integration starting at Tperi,2 is the reduction in the amplitudes of 
the eccentricity variations. 



both of the direct integrations shown in Figure 4, we have examined the 10 eccentricity-type mean- 
motion resonance variables at the 11:2 commensurability and confirmed that some of them alternate 
between circulation and libration. We have extended the direct integrations shown in Figure 4 and 
found that both are stable for at least 10^ yr, but we cannot rule out the possibility that the chaos 
would lead to instability on longer timescales. 

The orbital period P2 (and the other orbital parameters) of the outer planet of HD 12661 are not 
currently known to high precision, because the time span of the available observations is comparable 
to P2- To study the effects of varying P2, we have performed two sets of direct integrations with 
different initial P2, one starting at Tpgri,! and the other at rperi,2- The initial P2 were chosen such 
that (P2/Pi)/(ll/2) = 0.98, 0.99, 1.01, 1.02, and 1.03. The initial values of the other orbital 
parameters that can be obtained from the two-Kepler fit (Pi, i^i,2, ei^2, ^^1,2-, and rperi,i,2) were 
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HD 12661, sinj — 1, Jacobi coords 
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Fig. 5. — Same as Fig. 4, but for the sini = 1 HD 12661 system with an initial P2 such that 
P2/P1 = 0.99 X 11/2. This system shows regular secular evolution, with the direct integration 
results insensitive to the starting epoch. 



fixed at the values listed in Table 1, and the planetary masses and initial semimajor axes were 
derived assuming that sini = 1. We find that the cases with initial (P2/-Pi)/(ll/2) = 1.01 are also 
affected by the 11:2 mean-motion commensurability while those with initial {P2/Pi)/{ll/2) = 0.98, 
0.99, and 1.02 show regular secular evolution. Although the variations in aj for the cases with initial 
P2/P1 = 1.03 X 11/2 = 0.9997 x 17/3 indicate that these cases are probably affected by the 17:3 
commensurability, the variations in ej and the trajectories in the diagram of ei versus rui — VJ2 are 
qualitatively indistinguishable from those showing regular secular evolution (for at least 2 x 10^ yr). 

In Figures 5 and 6 wc show examples of the cases with regular secular evolution. Figure 5 shows 
the variations in aj and Cj for the sini = 1 HD 12661 system with P2/P1 = 0.99 x 11/2. The solid 
and dotted lines are from the direct integrations starting at Tperi,! and Tperi,2, respectively, and the 
dashed lines are from the octupole-level secular perturbation theory. The trajectories of this system 
in the phase-space diagrams of e\ versus w\ — W2 and 62 versus wi — W2 are the trajectories marked 
by the squares (which denote the current parameters of the HD 12661 system) in Figure 6. The 
direct integration results are not sensitive to the starting epoch. The amplitudes of eccentricity 
oscillations and the trajectories in the diagram of ei (or 62) versus vjx — W2 predicted by the 
octupole theory are in reasonably good agreement with those from direct integrations, even though 
a{= 0.323) is quite large. The main error of the octupole theory is in the period of eccentricity 
oscillations, with the predicted period (r^ 2.1 x 10^ yr) about 75% longer than that found in the 
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Fig. 6. — Same as Fig. 2, but for two-planet systems with the same masses, initial semimajor axes, 
and total angular momentum as the sinz = 1 HD 12661 system with initial P2/P1 = 0.99 x 11/2. 
The phase space of these systems with A(Ri 0.83) almost identical to Acrit(~ 0.82) is dominated 
by large libration islands about rui — W2 = 0° and 180°, and libration in the secular resonance is 
almost the only possibility. The trajectories through the squares are those of the HD 12661 system, 
with the squares showing the current parameters of the system. The initial conditions for the other 
trajectories are w\ — VJ2 = 0° and 62 = 0.01, 0.05, 0.10, 0.165, or 0.24, or wi — W2 = 180° and 
62 = 0.01, 0.05, 0.10 or 0.24, with ei being determined from the total angular momentum. All the 
direct integrations (except that for the HD 12661 system) start with the outer planet at apoapse 
and the inner planet at opposition. 

direct integrations (~ 1.2 x 10^ yr). 

Figures 5 and 6 show some interesting properties of the HD 12661 system that are shared 
by all of the cases studied above with different initial P2, including those affected by the close 
proximity to a mean-motion commensurability. The HD 12661 system is in a secular resonance 
with wi —W2 librating about 180°. This means that the lines of apsides of the two orbits are 
on average anti-aligned. However, the amplitude of libration of wi — W2 is large, and the orbital 
eccentricities of both planets exhibit large-amplitude variations. For the direct integrations shown 
in Figures 5 and 6, the amplitude of libration of zui — W2 is 56°, ei varies between 0.09 and 0.37, 
and 62 varies between 0.17 and 0.37. HD 12661 is the first extrasolar planetary system found to 
have vji — tU2 librating about 180° . The outer two planets of the v And system (Rivera k. Lissauer 
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2000; Lissauer & Rivera 2001; Chiang et al. 2001) and the two planets about GJ 876 (Laughhn & 
Chambers 2001; Lee &: Peale 2002) are also hkely in the secular resonance involving zui — W2, but 
in both of these cases wi — W2 libratcs about 0°. (In the case of GJ 876, the secular resonance is 
associated with the simultaneous librations of both lowest order mean-motion resonance variables 
at the 2:1 commensur ability.) 

The libration of vji — W2 and the large-amplitude variations of both ej in the HD 12661 system 
were anticipated by the analytic results derived in § 3.3. As we discussed in § 3.3, the octupole 
theory predicts that large libration islands and large-amplitude variations of both Sj are likely if 
A = L1/L2 ~ Acrit = 27^/(5 — 37^). The HD 12661 system has A pa 0.83, which is almost identical 
to Acrit ~ 0-82 for the dimensionless total angular momentum 7 ~ 0.96 of this system. In addition 
to the trajectories of the sini = 1 HD 12661 system with P2/P1 = 0.99 x 11/2, Figure 6 also 
shows the trajectories of systems with the same masses, initial semimajor axes, and total angular 
momentum. The direct-integration results (left panels) are in reasonably good agreement with the 
octupole results (right panels) and confirm that the phase space is dominated by large libration 
islands about zui — ■CU2 = 0° and 180°. The only discrepancy is a narrow region of circulation 
represented by a single trajectory in the direct-integration results. Thus libration of zui — VJ2 is 
almost the only possibility. We can see in Figure 6 that ei decreases from near its maximum 
possible value to near zero and 62 increases from near zero to near its maximum possible value at 
voi — W2 ^ 90°, and vice versa at wi — W2 ~ 270°, again in agreement with the analysis in § 3.3. 

Kiseleva-Eggleton et al. (2002) have recently reported that the sini = 1 HD 12661 system 
with the best-fit orbital parameters in Table 1 is regular based on the MEGNO (Mean Exponential 
Growth of Nearby Orbits) technique. However, they assumed that the best-fit orbital parameters 
are in astrocentric coordinates and showed results from only one (unspecified) starting epoch. We 
have repeated the direct integrations shown in Figure 4, but with the best-fit orbital parameters 
assumed to be in astrocentric coordinates. Because of the larger fiuctuations in the astrocentric 
02, the average ratio of the orbital periods is more sensitive to the starting epoch. The integration 
starting at Tperi,i turns out to be sufficiently far from the 11:2 commensurability that it shows 
regular secular evolution similar to those in Figure 5. The integration starting at rperi,2 is qualita- 
tively similar to the integration in Jacobi coordinates starting at the same epoch (dotted lines in 
Figure 4) and is thus cither regular or very weakly chaotic. These results again demonstrate the 
importance of interpreting the orbital parameters from multiple Kepler fits as orbital parameters 
in Jacobi coordinates. 

The best-fit orbital parameters of the HD 12661 planets in Table 1 were provided by D. A. Fis- 
cher (2002, private communication), and some of them differ slightly from those published later 
in Fischer et al. (2003). Adopting the parameters in Fischer et al. (2003) does not qualitatively 
change the results presented here, although the smaller w\ — W2 for the parameters in Fischer et al. 
(2003) does lead to slightly larger amplitude of libration of w\ — W2 and slightly larger variations 
in Cj. 
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After this paper was completed, we learned of contemporaneous work on the HD 12661 system 
by Gozdziewski (2003) and Gozdziewski & Maciejewski (2003). Gozdziewski (2003) has examined 
a wide region of the parameter space near the HD 12661 best-fit solution in Table 1, using as- 
trocentric coordinates and a combination of the MEGNO technique and direct integrations. His 
numerical results also show the proximity of the best-fit solution to the 11:2 commensurability and 
the robustness of the libration of wi — W2- Gozdziewski & Maciejewski (2003) have subsequently 
reported a best-fit solution to the data set published by Fischer et al. (2003) with P2 ~ 1660 days 
or P2/P1 ~ 6.3. This illustrates the uncertainties mentioned above in the orbital parameters of 
the outer planet due to the time span of the available observations being comparable to P2. Their 
best-fit solution also shows libration of tui — W2 about 180° and large-amplitude variations of both 
ej. 

Mayor et al. (2003) have reported the detection of two planets around the star HD 83443. 
There is considerable doubt about the existence of the outer planet in this system (Butler et al. 
2002), but it is interesting to note that the system as reported by Mayor et al. shares some of 
the properties of the HD 12661 system. The best-fit orbital periods of the HD 83443 planets arc 
nearly in the ratio 10:1, and direct integrations of this system (neglecting the effects of tides and 
general relativistic precession) show that its evolution is also affected by the close proximity to the 
mean-motion commensurability. The HD 83443 system also shows large-amplitude variations of 
both eccentricities, although its wi — VJ2 is circulating. The value of A(r^ 0.96) of the HD 83443 
system is not as close to its value of Acrit(~ 0.79) as in the HD 12661 system, and there is a larger 
region of circulation in the phase space of systems with the same masses, initial semimajor axes, 
and total angular momentum as the HD 83443 system. It should be noted that the inner planet of 
the HD 83443 system is sufficiently close to the star (P ~ 3 days) that tides and general relativistic 
precession can significantly affect the dynamics of the system (Wu & Goldreich 2002). 



4.3. Small Mutual Inclination 

To study the effects of a small mutual inclination angle imu between the orbits, we have repeated 
the direct integrations of the HD 168443 and HD 12661 systems shown in Figures 16 and 5 with 
initial imu = 1° or 2°. We assumed that the intersection of the orbital planes is initially along the 
line of sight, so that, initially, sinii = sini2 = 1 and the mutual inclination is the difference in 
the longitudes of ascending node referenced to the plane of the sky: = VL2 — ^i- The results 
are nearly identical to those with coplanar orbits, confirming the analysis in § 3.1 based on an 
expansion of the more general form of the octupole theory in powers of imu- 
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4.4. Stability 

The applicability of the octupole theory is limited by its use of an expansion in ri/r2 and orbit 
averaging. In particular, the orbit averaging eliminates the effects of mean-motion commensurabil- 
ities and the possible development of instabilities. The direct integrations reported in § 4.1 show 
that the coplanar HD 168443 system is stable for sin z > 0.4, and Marcy et al. (2001) have found 
that this system is stable for even smaller sini. The stability of the HD 168443 system is consis- 
tent with the empirical stability criteria for hierarchical triple systems derived most recently by 
Eggleton & Kiseleva (1995) and Mardhng & Aarseth (2001). For example, the sini = 1 HD 168443 
system has a = 01/02 = 0.102, while the stability criterion of Eggleton & Kiseleva requires ini- 
tial 02(1 — e2)/[ai(l + ei)] > 1.65 or a < 0.316 and that of Mardling & Aarseth requires initial 
02(1 — 62)701 > 3.17 or a < 0.252. The criterion of Eggleton & Kiseleva was tested over a wide 
range of mass ratios (including ones similar to those of the HD 168443 system) and is expected 
to be reliable to about 20%, while the criterion of Mardling & Aarseth does not have explicit de- 
pendence on mi/nio and may not be very accurate for mi /mo significantly different from unity. 
For a more detailed comparison with the above stability criteria, we have performed a set of direct 
integrations with different initial 02: 02 = 1.299-0.738 AU or a = 0.227-0.400. The other initial 
conditions are identical to those of the sini = 1 calculation shown as the solid lines in Figure 16. 
The evolutions of a arc shown in Figure 7. Except for the case with initial a = 0.333, which is 
apparently stabilized by the 5:1 mean- motion commensurability, there is a clear boundary between 
stable and unstable configurations at initial a ^ 0.30. We do not expect this stability boundary to 
move down significantly if we increase the integration length, since 10^ yr span many eccentricity 
oscillation cycles (if the orbits are stable) and, as we can see in Figure 7, instabilities usually man- 
ifest themselves on much shorter timescales. The stability criterion of Eggleton & Kiseleva is a few 
percent larger than the stability boundary determined here, consistent with the ~ 20% accuracy 
of the former. Even the stability criterion of Mardling &: Aarseth is only 16% smaller than the 
stability boundary determined here (see, however, the next paragraph). Finally, we note that some 
of the configurations in Figure 7 that are stable for at least 10^ yr appear to be chaotic (e.g., the 
case with initial a = 1/4). As in the case of the HD 12661 system, the chaos is due to the close 
proximity to a high-order mean- motion commensurability (8:1 for the case with initial a = 1/4). 

In § 4.2 we have considered the coplanar HD 12661 system with sini = 1 (i.e., minimum 
planetary masses) only. Unlike the HD 168443 system, which has a = 01/02 0.1 and is stable 
for any reasonable value of sini, the HD 12661 system has a ~ 0.32 (which is close to the stability 
boundary for the systems shown in Fig. 7 with the planetary masses of the sin i = 1 HD 168443 
system) and could be unstable for moderately small sini. To study the effects of sini, we have 
performed two sets of direct integrations with different sini. The initial values of Pj, Kj, ej, ujj, 
and Tpeii,j are the same as those of the HD 12661 system with P2/P1 = 0.99 x 11/2 shown in 
Figure 5, and the planetary masses and initial semimajor axes were derived for the assumed sini. 
The evolutions of a are shown in Figure 8 for sini = 0.1, 0.2,... ,1. The left and right panels are 
from the direct integrations starting at rperi,i and rperi,2; respectively. The dashed line in each 
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Fig. 7. — Evolutions of 01/02 = a for a set of direct integrations with different initial 02: 02 = 1.299- 
0.738 AU or a = 0.227-0.400. The other initial conditions are identical to those of the sini = 1 
HD 168443 calculation shown as the solid lines in Fig. lb. The dashed line in each panel indicates the 
initial a. Some of the cases that are stable for at least 10^ yr are affected by the close proximity to 
a high-order mean-motion commensurability (e.g., 8:1 for the case with initial a = 1/4) and appear 
to be chaotic. The case with initial a = 0.333 is apparently stabilized by the 5:1 commensurability. 



panel indicates the location of the 11:2 mean-motion commensurability. The influence of the 11:2 
commensurability increases with decreasing sini (which is expected since the width of a resonance 
increases with increasing planetary masses), but the cases with s'mi > 0.3 are stable (for at least 
2 X 10^ yr). All of the cases with sini < 0.3 are unstable. Are the stability criteria of Eggleton 
Sz Kiseleva (1995) and Mardling & Aarseth (2001) consistent with this transition from stable to 
unstable configurations at sinf ~ 0.3? For the configurations shown in Figure 8, which have the 
same initial Pi and P2, the initial a decreases slightly from 0.323 to 0.322 as sini decreases from 
1 to 0.1 (and mi and m2 increase; see eqs. [2] and [6]). For comparison, the stability boundary 
predicted by the criterion of Eggleton &: Kiseleva decreases from initial a = 0.444 to 0.348 as sini 
decreases from 1 to 0.1, with initial a = 0.397 at sinz = 0.3. This stability boundary is consistent 
with the transition at sini 0.3 if we simply move it down by 23%, which is compatible with 
the ~ 20% accuracy of the stability criterion of Eggleton & Kiseleva. On the other hand, the 
stability boundary predicted by the criterion of Mardling &; Aarseth changes only slightly from 
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Fig. 8.— Evolutions of ai/a2 = a for the HD 12661 system with P2/P1 = 0.99x 11/2 and sinz = 0.1, 
0.2, . . . , 1. The initial values of Pj, Kj, ej, ujj, and Tperij are the same as those of the system shown 
in Fig. 5, and the planetary masses and initial semimajor axes were derived for the assumed sini. 
The left and right panels are from direct integrations starting at Tpcri,i and Tpcri,25 respectively. 
The dashed line in each panel indicates the location of the 11:2 mean-motion commensur ability. 



initial a = 0.254 to 0.253 as sinz decreases from 1 to 0.1. Thus, this stability boundary is not 
consistent with the transition at sinz 0.3, unless we allow adjustments that change with sinz. 
This is not surprising since the stability criterion of Mardling & Aarseth does not have explicit 
dependence on mi/ mo and is not expected to be very accurate for mi /mo significantly different 
from unity. 



5. CONCLUSIONS 

We have investigated the dynamical evolution of coplanar hierarchical two-planet systems 
where the ratio of the orbital semimajor axes a = ai/a2 is small. Hierarchical two-planet systems 
are common among the known multiple planet systems and are likely to be common in the overall 
population of extrasolar planetary systems. We began by showing that the orbital parameters 
obtained from a two (or more generally multiple) Kepler fit to the radial velocity variations of a 
host star are best interpreted as Jacobi coordinates and that these coordinates should be used in any 
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analyses of hierarchical (and possibly other types of) planetary systems. We used the HD 168443 
system as an example to show that the use of astrocentric coordinates can introduce significant 
high-frequency variations in orbital elements that should be nearly constant on orbital timescales, 
leading to erroneous sensitivity of the evolution of the orbital elements on the starting epoch. 

An approximate theory that can be applied to hierarchical two-planet systems with a wide 
range of masses, orbital eccentricities, and inclinations is the octupole-level secular perturbation 
theory, which was developed for general hierarchical triple systems by Marchal (1990), Krymolowski 
k, Mazeh (1999), and Ford et al. (2000). The octupole theory is based on an expansion to order 
and orbit-averaging. We showed that the octupole approximation reduces the coplanar problem 
to one degree of freedom, with ei (or 62) and wi — W2 as the relevant phase-space variables. An 
analysis of the octupole equations yielded the following results: (1) the scaling properties of the 
equations imply that the amplitudes of eccentricity oscillations and the trajectory in the phase-space 
diagram of ei (or 62) versus wi — W2 should be independent of the inclination i of both orbits from 
the plane of the sky (as long as mi,m2 ^ mo) and that the period of eccentricity oscillations should 
be proportional to sini; (2) the fixed points in the phase-space diagram must be at zui — W2 = 0° or 
180°; and (3) if the ratio of the maximum orbital angular momenta, A = L1/L2 ~ (mi/m2)a^/^, for 
given semimajor axes is approximately equal to Acrit = 27^/(5 — 37^), where 7 = {G1+G2) / {L1+L2) 
is the dimensionless total angular momentum, then libration of wi — W2 is almost certain, with 
possibly large amplitude variations of both eccentricities. 

We used both the octupole-level secular perturbation theory and direct numerical orbit in- 
tegrations to study the dynamical evolution of the HD 168443 and HD 12661 systems and their 
variants. The HD 168443 system is not in a secular resonance and its wi — W2 circulates. For the 
family of systems with the same masses, initial semimajor axes, and total angular momentum as 
the smi = 1 HD 168443 system, there are two relatively small libration islands in the phase-space 
diagram of ei (or 62) versus wi — W2-, one about a fixed point at vji — VJ2 = 0° and the other 
about a fixed point at wi — W2 = 180°, and the trajectory of the HD 168443 system is far from 
both. In all cases, the octupole results are in excellent agreement with the direct-integration results. 
Direct integrations of the HD 168443 system with sini = 0.4 also confirmed the analytic octupole 
results on the effects of sini. To demonstrate how the octupole theory can be used to explore the 
parameter space rapidly, we used the octupole theory to determine the number and positions of 
fixed points as a function of the total angular momentum for systems with the same masses and 
semimajor axes as the sini = 1 HD 168443 system. For sufficiently low total angular momentum, 
the existence of fixed point (s) in addition to those appropriate for the total angular momentum of 
the HD 168443 system leads to systems with unusual behaviors such as e\ increasing to unity in 
a finite time or libration of wi — VJ2 about 0° with e\ very close to 1. These results were again 
confirmed by direct integrations. Direct integrations of systems similar to the sini = 1 HD 168443 
system, but with different initial a,2, showed that these systems with relatively massive planets 
[(mi -|- m2)/mo ~ 0.024] are unstable if the initial a. > 0.30. 

The HD 12661 system is the first extrasolar planetary system found to have tui — W2 librating 
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about 180°. The secular resonance means that the lines of apsides of the two orbits are on average 
ant i- aligned, although the amphtude of hbration of w\ — is large. The libration of vox — W2 
and the large-amplitude variations of both orbital eccentricities in the HD 12661 system, which 
has A « 0.83 and Acrit ~ 0.82, are consistent with the analytic results on systems with A « 
Acrit- Direct integrations showed that the evolution of the sini = 1 HD 12661 system with the 
best-fit orbital parameters (-P2/-F1 = 0.9975 x 11/2) is affected by the close proximity to the 
11:2 mean-motion commensurability, but that small changes in the orbital period of the outer 
planet within the uncertainty can result in configurations that are not affected by mean-motion 
commcnsurabilities. For the family of systems with the same masses, initial semimajor axes, and 
total angular momentum as the sini = 1 HD 12661 system with P2/P1 = 0.99 x 11/2, which are not 
affected by mean-motion commensur abilities, the octupole results are in reasonably good agreement 
with the direct-integration results, even though a{'^ 0.32) is quite large. The phase space of this 
family of systems with A ~ Acrit is dominated by large libration islands about wi — vj2 = 0° and 
180°, which again confirmed the analytic results on systems with A ~ Acrit j and libration in the 
secular resonance is almost the only possibility. Finally, we showed that the HD 12661 system is 
unstable if sini < 0.3 [or (mi -|- 1112) /rriQ > 0.012]. 

The comparisons with direct numerical orbit integrations established that the octupole-level 
secular perturbation theory is highly accurate for coplanar hierarchical two-planet systems with 
a < 0.1 and reasonably accurate even for systems with a as large as 1/3, provided that a is 
not too close to a significant mean-motion commensurability or above the stability boundary. We 
have focused our study in this paper on coplanar systems in which the star and the planets can 
be treated as point masses interacting via Newtonian gravity, where a small mutual inclination 
between the orbits is shown to have little effect on the secular evolution. The more general form 
of the octupole theory is applicable to systems with arbitrary mutual inclinations. If the inner or 
both planets are very close to the star, tidal interactions and general relativistic precession can be 
important (see, e.g., Wu Sz Goldreich 2002). Nagasawa et al. (2003) have also discussed the effects 
of additional orbital precessions due to interactions with the protoplanetary disk during the epoch 
of disk depletion. It is straightforward to modify the octupole theory to include any of these effects 
(see, e.g., Blaes et al. 2002 for relativistic precession). 

We thank Debra Fischer for furnishing the orbital parameters of the HD 12661 system before 
publication. We also thank O. Blaes, G. W. Marcy, N. Murray, E. J. Rivera, and Y. Wu for infor- 
mative discussions and an anonymous referee for comments that greatly improved the presentation. 
This research was supported in part by NASA grants NAG5 11666 and NAG5 12087. 
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